Setup

library(magrittr)
library(dplyr)
library(ggplot2)
library(cowplot)
library(ggrepel)
library(ggpubr)
library(scHelper)
library(biomaRt)
library(pheatmap)
library(SummarizedExperiment)
library(qs)
system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(DEP)

tt <- Sys.time()

Load data

tmp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")

dat.clean.orig <- tmp$dat.clean
df.clean <- tmp$expDesign.clean
universe <- tmp$universe

Metadata

# Metadata
metadata <- read.delim("/work/proteomics/02_data/Metadata_Cohort_100723.csv", sep = ";", dec = ",", header = T) %>% 
  mutate(intern = gsub("K", "", .$K.no.),
         Diagnosis = gsub("Unknown", "MSA", .$Diagnosis)) # 1418 should be MSA

dat.clean.orig@colData$subtype <- metadata$Subtype[match(dat.clean.orig@colData$intern, metadata$intern)]

REGION-SPECIFIC PROTEINS WITH IMPUTATION

Clean data

imp <- dat.clean.orig@elementMetadata@listData %>% 
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

dat.tmp <- dat.clean.orig %>% 
  .[rownames(.) %in% rownames(imp)]

Calculate DEPs

# dat.tmp <- dat.clean
dat.tmp@elementMetadata$ID <- dat.tmp@elementMetadata$CER_ID
dat.tmp@elementMetadata$name <- dat.tmp@elementMetadata$CER_name
all_contrasts <- test_diff(dat.tmp, type = "all")
## Tested contrasts: CER_vs_FC, CER_vs_MED, CER_vs_STR, CER_vs_SN, FC_vs_MED, FC_vs_STR, FC_vs_SN, MED_vs_STR, MED_vs_SN, STR_vs_SN
dep_all <- add_rejections(all_contrasts, alpha = 0.05, lfc = 2)

Number of significant proteins

comb <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

comb %>% 
  lapply(\(x) dep_all@elementMetadata %>% 
           as.data.frame() %>% 
           .[, paste(x, "significant", sep = "_")]) %>% 
  setNames(comb) %>% 
  sapply(sum) %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(line = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "", y = "No. significant proteins", fill = "Comparison") +
  scale_fill_manual(values = ggsci::pal_cosmic()(10)) +
  geom_text(aes(y = value + 8, label = value))

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, 45)) +
  ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Correlation

plot_cor(dep_all, significant = TRUE, lower = 0.6, upper = 1, pal = "YlOrRd")

Clustering of DEPs per sample

ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE, 
             k = 3, col_limit = 3, show_row_names = FALSE,
             indicate = c("condition"))

protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()

rowIds <- ComplexHeatmap::row_order(ht) %>% 
  lapply(\(x) protIds[x])

go.res <- rowIds %>% 
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
## 
## 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

# No significant terms, instead we ouput top hits
go.res %>% 
  lapply(getElement, "result") %>%
  lapply(dplyr::slice, seq(10))
## $`3`
##                    ID
## GO:0090263 GO:0090263
## GO:0060828 GO:0060828
## GO:0030111 GO:0030111
## GO:0003018 GO:0003018
## GO:0071260 GO:0071260
## GO:0071496 GO:0071496
## GO:0060070 GO:0060070
## GO:0090288 GO:0090288
## GO:0051882 GO:0051882
## GO:0006972 GO:0006972
##                                                                   Description
## GO:0090263             positive regulation of canonical Wnt signaling pathway
## GO:0060828                      regulation of canonical Wnt signaling pathway
## GO:0030111                                regulation of Wnt signaling pathway
## GO:0003018                             vascular process in circulatory system
## GO:0071260                           cellular response to mechanical stimulus
## GO:0071496                             cellular response to external stimulus
## GO:0060070                                    canonical Wnt signaling pathway
## GO:0090288 negative regulation of cellular response to growth factor stimulus
## GO:0051882                                       mitochondrial depolarization
## GO:0006972                                              hyperosmotic response
##            GeneRatio  BgRatio RichFactor FoldEnrichment   zScore      pvalue
## GO:0090263     6/191  46/6858 0.13043478       4.683360 4.242328 0.001597978
## GO:0060828     9/191 100/6858 0.09000000       3.231518 3.804608 0.001785780
## GO:0030111    10/191 123/6858 0.08130081       2.919167 3.635090 0.002172139
## GO:0003018    11/191 146/6858 0.07534247       2.705228 3.524939 0.002427273
## GO:0071260     5/191  35/6858 0.14285714       5.129394 4.145255 0.002616427
## GO:0071496     5/191  35/6858 0.14285714       5.129394 4.145255 0.002616427
## GO:0060070     9/191 107/6858 0.08411215       3.020111 3.564524 0.002849182
## GO:0090288     5/191  36/6858 0.13888889       4.986911 4.059294 0.002970923
## GO:0051882     3/191  12/6858 0.25000000       8.976440 4.680581 0.003885627
## GO:0006972     4/191  25/6858 0.16000000       5.744921 4.022646 0.004661548
##             p.adjust    qvalue
## GO:0090263 0.8816214 0.8816214
## GO:0060828 0.8816214 0.8816214
## GO:0030111 0.8816214 0.8816214
## GO:0003018 0.8816214 0.8816214
## GO:0071260 0.8816214 0.8816214
## GO:0071496 0.8816214 0.8816214
## GO:0060070 0.8816214 0.8816214
## GO:0090288 0.8816214 0.8816214
## GO:0051882 0.9169262 0.9169262
## GO:0006972 0.9169262 0.9169262
##                                                                         geneID
## GO:0090263                                 GPC5/PPM1B/PPM1A/USP8/USP47/TBL1XR1
## GO:0060828                 GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/USP47/HECW1/TBL1XR1
## GO:0030111           GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/RACK1/USP47/HECW1/TBL1XR1
## GO:0003018 SLC4A3/SLC44A1/SOD3/PRKG1/SLC38A2/FAAH/MANF/GCLC/AKAP12/SLC2A1/TJP2
## GO:0071260                                     FAS/SLC38A2/TMEM87A/GCLC/SLC2A1
## GO:0071496                                     FAS/SLC38A2/TMEM87A/GCLC/SLC2A1
## GO:0060070                 GPC5/DKK3/PPM1B/PPM1A/RBX1/USP8/USP47/HECW1/TBL1XR1
## GO:0090288                                           HRG/CADM4/PPM1A/EPN2/CASK
## GO:0051882                                                      ALB/GCLC/RACK1
## GO:0006972                                        RCSD1/SLC25A23/ICOSLG/SLC2A1
##            Count
## GO:0090263     6
## GO:0060828     9
## GO:0030111    10
## GO:0003018    11
## GO:0071260     5
## GO:0071496     5
## GO:0060070     9
## GO:0090288     5
## GO:0051882     3
## GO:0006972     4
## 
## $`2`
##                    ID                                          Description
## GO:0006897 GO:0006897                                          endocytosis
## GO:0072350 GO:0072350                 tricarboxylic acid metabolic process
## GO:0006886 GO:0006886                      intracellular protein transport
## GO:0043523 GO:0043523               regulation of neuron apoptotic process
## GO:0021955 GO:0021955           central nervous system neuron axonogenesis
## GO:0043547 GO:0043547               positive regulation of GTPase activity
## GO:0001845 GO:0001845                               phagolysosome assembly
## GO:0048169 GO:0048169 regulation of long-term neuronal synaptic plasticity
## GO:0043114 GO:0043114                  regulation of vascular permeability
## GO:0060627 GO:0060627             regulation of vesicle-mediated transport
##            GeneRatio  BgRatio RichFactor FoldEnrichment   zScore       pvalue
## GO:0006897    26/180 444/6858 0.05855856       2.231081 4.403457 8.400259e-05
## GO:0072350     4/180  12/6858 0.33333333      12.700000 6.659447 1.926603e-04
## GO:0006886    24/180 437/6858 0.05491991       2.092449 3.874543 4.224967e-04
## GO:0043523    10/180 125/6858 0.08000000       3.048000 3.793681 1.576887e-03
## GO:0021955     4/180  20/6858 0.20000000       7.620000 4.867305 1.599618e-03
## GO:0043547     7/180  70/6858 0.10000000       3.810000 3.879400 2.270127e-03
## GO:0001845     3/180  11/6858 0.27272727      10.390909 5.117218 2.511748e-03
## GO:0048169     4/180  23/6858 0.17391304       6.626087 4.436916 2.749374e-03
## GO:0043114     4/180  24/6858 0.16666667       6.350000 4.310246 3.232495e-03
## GO:0060627    19/180 366/6858 0.05191257       1.977869 3.156545 3.275084e-03
##             p.adjust    qvalue
## GO:0006897 0.2171467 0.2154003
## GO:0072350 0.2490135 0.2470108
## GO:0006886 0.3640513 0.3611235
## GO:0043523 0.7573412 0.7512504
## GO:0021955 0.7573412 0.7512504
## GO:0043547 0.7573412 0.7512504
## GO:0001845 0.7573412 0.7512504
## GO:0048169 0.7573412 0.7512504
## GO:0043114 0.7573412 0.7512504
## GO:0060627 0.7573412 0.7512504
##                                                                                                                                                                   geneID
## GO:0006897 SH3GL3/AP1S1/AP2S1/ACTG1/TSC2/SGIP1/YES1/ARL8B/EPS15L1/CSNK1G1/REPS1/SH3GL2/PPP3R1/PPP3CA/LYST/CORO1A/PIK3C3/SNX12/BCR/APP/RALB/RAB1A/HRAS/VAC14/PRKCG/VPS33B
## GO:0072350                                                                                                                                        IDH3A/IDH3B/IREB2/IDH2
## GO:0006886            STXBP1/AP1S1/AP2S1/TSC2/TOMM22/TOMM70/PPP3R1/PPP3CA/RAB24/IPO13/NF1/SEC24B/COPZ1/PIK3C3/XPO7/BCR/APBA1/NDEL1/TIMM13/RAB1A/STX7/USP9X/EXOC6B/VPS33B
## GO:0043523                                                                                                      STXBP1/PHB1/DNAJC5/CORO1A/NF1/AARS1/SARM1/APP/HRAS/PRKCG
## GO:0021955                                                                                                                                      ARHGAP35/DCC/KIFBP/NDEL1
## GO:0043547                                                                                                                           RGS6/RALGAPB/NF1/BCR/NDEL1/HRAS/ABR
## GO:0001845                                                                                                                                           ARL8B/CORO1A/VPS33B
## GO:0048169                                                                                                                                             FXR1/NF1/APP/HRAS
## GO:0043114                                                                                                                                      ARHGAP35/YES1/SH3GL2/BCR
## GO:0060627                                           STXBP1/SH3GL3/AP2S1/ACTG1/TSC2/SGIP1/CHMP3/PPP3R1/PPP3CA/DNAJC5/RAB3C/PRKCB/CORO1A/SNX12/BCR/APP/VAC14/PRKCG/TBC1D4
##            Count
## GO:0006897    26
## GO:0072350     4
## GO:0006886    24
## GO:0043523    10
## GO:0021955     4
## GO:0043547     7
## GO:0001845     3
## GO:0048169     4
## GO:0043114     4
## GO:0060627    19
## 
## $`1`
##                    ID                                            Description
## GO:0006397 GO:0006397                                        mRNA processing
## GO:0016071 GO:0016071                                 mRNA metabolic process
## GO:0042776 GO:0042776 proton motive force-driven mitochondrial ATP synthesis
## GO:0051028 GO:0051028                                         mRNA transport
## GO:0015986 GO:0015986               proton motive force-driven ATP synthesis
## GO:0007158 GO:0007158                              neuron cell-cell adhesion
## GO:0006403 GO:0006403                                       RNA localization
## GO:0008380 GO:0008380                                           RNA splicing
## GO:0050657 GO:0050657                                 nucleic acid transport
## GO:0050658 GO:0050658                                          RNA transport
##            GeneRatio  BgRatio RichFactor FoldEnrichment   zScore       pvalue
## GO:0006397    16/128 318/6858 0.05031447       2.695755 4.270240 0.0002583117
## GO:0016071    19/128 445/6858 0.04269663       2.287605 3.873439 0.0005451515
## GO:0042776     6/128  60/6858 0.10000000       5.357813 4.675393 0.0008247422
## GO:0051028     7/128  84/6858 0.08333333       4.464844 4.406210 0.0009203638
## GO:0015986     6/128  65/6858 0.09230769       4.945673 4.407695 0.0012612538
## GO:0007158     3/128  13/6858 0.23076923      12.364183 5.655717 0.0015840344
## GO:0006403     8/128 131/6858 0.06106870       3.271947 3.620659 0.0030021369
## GO:0008380    13/128 292/6858 0.04452055       2.385327 3.336245 0.0030064054
## GO:0050657     7/128 105/6858 0.06666667       3.571875 3.662357 0.0033499474
## GO:0050658     7/128 105/6858 0.06666667       3.571875 3.662357 0.0033499474
##             p.adjust    qvalue
## GO:0006397 0.3686057 0.3686057
## GO:0016071 0.3686057 0.3686057
## GO:0042776 0.3686057 0.3686057
## GO:0051028 0.3686057 0.3686057
## GO:0015986 0.4041057 0.4041057
## GO:0007158 0.4229372 0.4229372
## GO:0006403 0.4482364 0.4482364
## GO:0008380 0.4482364 0.4482364
## GO:0050657 0.4482364 0.4482364
## GO:0050658 0.4482364 0.4482364
##                                                                                                                            geneID
## GO:0006397                     SNRNP40/DDX42/ADAR/CSDC2/PNN/RBM25/ARVCF/SAFB/TRA2B/CSTF3/HNRNPF/SNRNP70/DDX39B/THOC6/CD2BP2/SRRM1
## GO:0016071 SNRNP40/DDX42/CNOT10/ADAR/CSDC2/PNN/RBM25/ARVCF/SAFB/TRA2B/CSTF3/CNOT9/HNRNPF/GIGYF2/SNRNP70/DDX39B/THOC6/CD2BP2/SRRM1
## GO:0042776                                                                               SDHD/NDUFB1/NDUFA11/NDUFB9/NDUFA8/ATP5MF
## GO:0051028                                                                           NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
## GO:0015986                                                                               SDHD/NDUFB1/NDUFA11/NDUFB9/NDUFA8/ATP5MF
## GO:0007158                                                                                                      CNTN4/NLGN3/ASTN2
## GO:0006403                                                                     NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2/STAU2
## GO:0008380                                      SNRNP40/DDX42/PNN/RBM25/ARVCF/TRA2B/HNRNPF/SNRNP70/DDX39B/RTCB/THOC6/CD2BP2/SRRM1
## GO:0050657                                                                           NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
## GO:0050658                                                                           NUP58/SEH1L/DDX39B/CHTOP/NUP153/THOC6/RANBP2
##            Count
## GO:0006397    16
## GO:0016071    19
## GO:0042776     6
## GO:0051028     7
## GO:0015986     6
## GO:0007158     3
## GO:0006403     8
## GO:0008380    13
## GO:0050657     7
## GO:0050658     7

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = TRUE, 
             k = 3, col_limit = 5, show_row_names = FALSE)

Volcano plots

combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1
## Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V2
## Warning: ggrepel: 58 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V3
## Warning: ggrepel: 161 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V4
## Warning: ggrepel: 137 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5

## 
## $V6

## 
## $V7

## 
## $V8
## Warning: ggrepel: 48 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V9
## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V10

Extract results

res <- dep_all@elementMetadata@listData %>% 
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x)))

# res2 <- get_results(dep_all)

res.split <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>% filter(!!sym(paste0(cc, "_p.adj")) <= 0.05, abs(!!sym(paste0(cc, "_diff"))) >= 1.5)), .)

GO on DEPs

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: KPNA2,APLF,HELQ,KAT5,SMARCAD1,MEAF6
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: SLC25A36,SLC5A1,SPIDR,STAT6,IL4,RRP7BP
## --> return NULL...
res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: NEIL3,H1-8,SLC25A33,KDM1A,EIF6,CLCF1
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: APTX,H1-7,APLF,TERF2IP,SESN2,NDFIP1
## --> return NULL...
res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Get markers

all.markers <- res.split %>% 
  names() %>% 
  lapply(\(comp) {
    up <- strsplit(comp, "_vs_") %>% sget(1)
    down <- strsplit(comp, "_vs_") %>% sget(2)
    diff <- paste0(comp, "_diff")
    
    up.df <- res.split[[comp]] %>% 
      filter(!!sym(diff) > 0) %>% 
      mutate(area = up) %>% 
      dplyr::select(name, !!sym(diff), area)
    
    down.df <- res.split[[comp]] %>% 
      filter(!!sym(diff) < 0) %>% 
      mutate(area = down) %>% 
      dplyr::select(name, !!sym(diff), area)
    
    out <- rbind(up.df, down.df) %>% 
      dplyr::rename(l2fc = !!sym(diff)) %>% 
      mutate(l2fc = abs(l2fc))
    
    return(out)
  }) %>% 
  bind_rows()

markers.sum <- all.markers %>% 
  dplyr::group_by(area, name) %>% 
  summarize(n = n(),
            mean_l2fc = mean(l2fc),
            max_l2fc = max(l2fc)) %>%
  data.frame()
## `summarise()` has grouped output by 'area'. You can override using the
## `.groups` argument.
markers.split <- markers.sum %>%  
  arrange(-n, -mean_l2fc) %$% 
  split(., area)

Plot top markers

markers.split %>% 
  lapply(dplyr::slice, seq(10)) %>% 
  lapply(pull, name) %>%
  Map(\(genes, name) plot_single(dep = dep_all, proteins = genes, type = "centered") + guides(col = "none") + labs(title = name), genes = ., name = names(.))
## $CER

## 
## $FC

## 
## $MED

## 
## $SN

## 
## $STR

Onthology on top markers

markers.enrich <- markers.split %>% 
  lapply(dplyr::slice, seq(50)) %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

markers.enrich %>% 
  lapply(head, 20) %>% 
  lapply(knitr::kable)
## $CER
## 
## 
## |           |ID         |Description                                                                          |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|    zScore|   pvalue|  p.adjust|    qvalue|geneID                                                                                                                           | Count|
## |:----------|:----------|:------------------------------------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|--------:|---------:|---------:|:--------------------------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0006397 |GO:0006397 |mRNA processing                                                                      |19/49     |318/6858 |  0.0597484|       8.362341| 11.404180| 0.00e+00| 0.0000000| 0.0000000|WDR33/ALKBH5/SRSF2/SNU13/SNRPD3/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1         |    19|
## |GO:0016071 |GO:0016071 |mRNA metabolic process                                                               |21/49     |445/6858 |  0.0471910|       6.604815| 10.371323| 0.00e+00| 0.0000000| 0.0000000|WDR33/ALKBH5/SRSF2/SNU13/FUS/FTO/SNRPD3/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1 |    21|
## |GO:0000377 |GO:0000377 |RNA splicing, via transesterification reactions with bulged adenosine as nucleophile |16/49     |215/6858 |  0.0744186|      10.415567| 11.898906| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1                            |    16|
## |GO:0000398 |GO:0000398 |mRNA splicing, via spliceosome                                                       |16/49     |215/6858 |  0.0744186|      10.415567| 11.898906| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1                            |    16|
## |GO:0000375 |GO:0000375 |RNA splicing, via transesterification reactions                                      |16/49     |219/6858 |  0.0730594|      10.225328| 11.769988| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1                            |    16|
## |GO:1903241 |GO:1903241 |U2-type prespliceosome assembly                                                      |8/49      |22/6858  |  0.3636364|      50.894249| 19.883107| 0.00e+00| 0.0000000| 0.0000000|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF                                                                               |     8|
## |GO:0008380 |GO:0008380 |RNA splicing                                                                         |17/49     |292/6858 |  0.0582192|       8.148309| 10.589312| 0.00e+00| 0.0000000| 0.0000000|SRSF2/SNU13/FUS/SNRPD3/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF7/SRSF3/SNRPA1/SNRPE/SNRNP70/SNRPF/USP39/KHDRBS1                        |    17|
## |GO:0000245 |GO:0000245 |spliceosomal complex assembly                                                        |9/49      |60/6858  |  0.1500000|      20.993878| 13.194882| 0.00e+00| 0.0000000| 0.0000000|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39                                                                         |     9|
## |GO:0036260 |GO:0036260 |RNA capping                                                                          |5/49      |12/6858  |  0.4166667|      58.316327| 16.856755| 0.00e+00| 0.0000010| 0.0000009|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF                                                                                                  |     5|
## |GO:0022618 |GO:0022618 |protein-RNA complex assembly                                                         |9/49      |134/6858 |  0.0671642|       9.400244|  8.330170| 3.00e-07| 0.0000273| 0.0000249|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39                                                                         |     9|
## |GO:0000966 |GO:0000966 |RNA 5'-end processing                                                                |5/49      |23/6858  |  0.2173913|      30.425909| 11.990796| 5.00e-07| 0.0000341| 0.0000311|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF                                                                                                  |     5|
## |GO:0071826 |GO:0071826 |protein-RNA complex organization                                                     |9/49      |141/6858 |  0.0638298|       8.933565|  8.074464| 5.00e-07| 0.0000351| 0.0000321|SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39                                                                         |     9|
## |GO:0000387 |GO:0000387 |spliceosomal snRNP assembly                                                          |5/49      |28/6858  |  0.1785714|      24.992711| 10.791242| 1.30e-06| 0.0000820| 0.0000748|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF                                                                                                  |     5|
## |GO:0001510 |GO:0001510 |RNA methylation                                                                      |5/49      |31/6858  |  0.1612903|      22.574062| 10.212246| 2.20e-06| 0.0001295| 0.0001182|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF                                                                                                  |     5|
## |GO:0045934 |GO:0045934 |negative regulation of nucleobase-containing compound metabolic process              |14/49     |474/6858 |  0.0295359|       4.133816|  5.998458| 3.30e-06| 0.0001808| 0.0001650|H1-10/H1-4/HMGB1/ZNF512B/HDGF/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/SRSF7/CCAR1/KHDRBS1                                               |    14|
## |GO:0022613 |GO:0022613 |ribonucleoprotein complex biogenesis                                                 |10/49     |238/6858 |  0.0420168|       5.880638|  6.500691| 5.20e-06| 0.0002528| 0.0002307|SNU13/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1/SNRPE/SNRPF/USP39                                                                   |    10|
## |GO:0051253 |GO:0051253 |negative regulation of RNA metabolic process                                         |13/49     |423/6858 |  0.0307329|       4.301346|  5.945801| 5.30e-06| 0.0002528| 0.0002307|H1-4/HMGB1/ZNF512B/HDGF/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/SRSF7/CCAR1/KHDRBS1                                                     |    13|
## |GO:0043484 |GO:0043484 |regulation of RNA splicing                                                           |7/49      |117/6858 |  0.0598291|       8.373626|  6.823942| 1.67e-05| 0.0007526| 0.0006869|SRSF2/FUS/SF3B3/SRSF7/SRSF3/SNRNP70/KHDRBS1                                                                                      |     7|
## |GO:0043414 |GO:0043414 |macromolecule methylation                                                            |5/49      |49/6858  |  0.1020408|      14.281549|  7.914593| 2.28e-05| 0.0009724| 0.0008874|SNRPD3/SNRPD2/SNRPB/SNRPE/SNRPF                                                                                                  |     5|
## |GO:1903311 |GO:1903311 |regulation of mRNA metabolic process                                                 |8/49      |178/6858 |  0.0449438|       6.290300|  6.066323| 3.20e-05| 0.0012974| 0.0011841|ALKBH5/SRSF2/FUS/FTO/SRSF7/SRSF3/SNRNP70/KHDRBS1                                                                                 |     8|
## 
## $FC
## 
## 
## |           |ID         |Description                                     |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|   zScore|    pvalue|  p.adjust|    qvalue|geneID                                                                            | Count|
## |:----------|:----------|:-----------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:---------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission                  |13/47     |495/6858 |  0.0262626|       3.832108| 5.433655| 0.0000178| 0.0076320| 0.0069256|EPHA4/CAMKV/SLC17A7/CHRM1/CAMK2A/NPY/GABRA1/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A |    13|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling            |13/47     |495/6858 |  0.0262626|       3.832108| 5.433655| 0.0000178| 0.0076320| 0.0069256|EPHA4/CAMKV/SLC17A7/CHRM1/CAMK2A/NPY/GABRA1/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A |    13|
## |GO:0021954 |GO:0021954 |central nervous system neuron development       |5/47      |51/6858  |  0.0980392|      14.305382| 7.922187| 0.0000226| 0.0076320| 0.0069256|PLXNA1/EPHA4/GABRB1/NPY/SLC4A10                                                   |     5|
## |GO:0050808 |GO:0050808 |synapse organization                            |11/47     |394/6858 |  0.0279188|       4.073766| 5.220100| 0.0000521| 0.0132177| 0.0119943|SYNPO/PLXNA1/EPHA4/ACTN1/ICAM5/CAMKV/LINGO2/GABRA1/AMIGO1/TNC/HOMER1              |    11|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission    |10/47     |346/6858 |  0.0289017|       4.217193| 5.101156| 0.0000911| 0.0157871| 0.0143259|EPHA4/CAMKV/CHRM1/CAMK2A/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A                    |    10|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling          |10/47     |347/6858 |  0.0288184|       4.205040| 5.089615| 0.0000933| 0.0157871| 0.0143259|EPHA4/CAMKV/CHRM1/CAMK2A/VGF/PRKAR2B/SV2B/SLC4A10/HOMER1/STX1A                    |    10|
## |GO:0060078 |GO:0060078 |regulation of postsynaptic membrane potential   |5/47      |80/6858  |  0.0625000|       9.119681| 6.067974| 0.0001990| 0.0283848| 0.0257575|GABRB1/CHRM1/GABRA1/RGS7BP/STX1A                                                  |     5|
## |GO:0060560 |GO:0060560 |developmental growth involved in morphogenesis  |6/47      |129/6858 |  0.0465116|       6.786739| 5.511427| 0.0002249| 0.0283848| 0.0257575|PLXNA1/ISLR2/CPNE5/RASAL1/OLFM1/TNC                                               |     6|
## |GO:0007186 |GO:0007186 |G protein-coupled receptor signaling pathway    |9/47      |318/6858 |  0.0283019|       4.129667| 4.747154| 0.0002531| 0.0283848| 0.0257575|NECAB2/CHRM1/CAMK2A/NPY/RGS7BP/CHGA/PRKAR2B/GNG2/HOMER1                           |     9|
## |GO:0021953 |GO:0021953 |central nervous system neuron differentiation   |5/47      |86/6858  |  0.0581395|       8.483424| 5.800984| 0.0002797| 0.0283848| 0.0257575|PLXNA1/EPHA4/GABRB1/NPY/SLC4A10                                                   |     5|
## |GO:0006821 |GO:0006821 |chloride transport                              |4/47      |56/6858  |  0.0714286|      10.422492| 5.881012| 0.0005480| 0.0459986| 0.0417410|GABRB1/SLC17A7/GABRA1/SLC4A10                                                     |     4|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport |5/47      |101/6858 |  0.0495050|       7.223510| 5.233951| 0.0005891| 0.0459986| 0.0417410|SCN3B/CHRM1/CAMK2A/AMIGO1/HOMER1                                                  |     5|
## |GO:1990138 |GO:1990138 |neuron projection extension                     |5/47      |101/6858 |  0.0495050|       7.223510| 5.233951| 0.0005891| 0.0459986| 0.0417410|PLXNA1/ISLR2/CPNE5/RASAL1/OLFM1                                                   |     5|
## 
## $MED
## 
## 
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
## 
## $SN
## 
## 
## |           |ID         |Description                                 |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|   zScore|    pvalue|  p.adjust|    qvalue|geneID                                                         | Count|
## |:----------|:----------|:-------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:--------------------------------------------------------------|-----:|
## |GO:0042220 |GO:0042220 |response to cocaine                         |5/46      |38/6858  |  0.1315789|      19.616705| 9.456079| 0.0000046| 0.0045342| 0.0040917|SLC6A4/PPP1R1B/GAD2/HTR1B/CNR1                                 |     5|
## |GO:0043279 |GO:0043279 |response to alkaloid                        |5/46      |65/6858  |  0.0769231|      11.468227| 6.967993| 0.0000665| 0.0183148| 0.0165274|SLC6A4/PPP1R1B/GAD2/HTR1B/CNR1                                 |     5|
## |GO:0006837 |GO:0006837 |serotonin transport                         |3/46      |13/6858  |  0.2307692|      34.404682| 9.906060| 0.0000771| 0.0183148| 0.0165274|SLC6A4/HTR1B/CNR1                                              |     3|
## |GO:0007631 |GO:0007631 |feeding behavior                            |4/46      |37/6858  |  0.1081081|      16.117509| 7.576450| 0.0000995| 0.0183148| 0.0165274|PMCH/HTR1B/CNR1/CARTPT                                         |     4|
## |GO:0003013 |GO:0003013 |circulatory system process                  |9/46      |292/6858 |  0.0308219|       4.595146| 5.159008| 0.0001115| 0.0183148| 0.0165274|SLC2A13/SLC6A4/SCN3B/HTR1B/ATP8A1/CNR1/KCNMA1/CARTPT/SCN2B     |     9|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport      |8/46      |227/6858 |  0.0352423|       5.254166| 5.356073| 0.0001120| 0.0183148| 0.0165274|SLC6A4/TESC/DPP6/SCN3B/HTR1B/CNR1/HPCA/SCN2B                   |     8|
## |GO:0015844 |GO:0015844 |monoamine transport                         |4/46      |45/6858  |  0.0888889|      13.252174| 6.775774| 0.0002159| 0.0289610| 0.0261347|SLC6A4/HTR1B/CNR1/CARTPT                                       |     4|
## |GO:0007617 |GO:0007617 |mating behavior                             |3/46      |19/6858  |  0.1578947|      23.540046| 8.084329| 0.0002539| 0.0289610| 0.0261347|SLC6A4/PPP1R1B/CNR1                                            |     3|
## |GO:0023061 |GO:0023061 |signal release                              |8/46      |260/6858 |  0.0307692|       4.587291| 4.845692| 0.0002857| 0.0289610| 0.0261347|SLC6A4/VGF/WLS/SYT5/SLC32A1/HTR1B/CNR1/CARTPT                  |     8|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission              |11/46     |495/6858 |  0.0222222|       3.313044| 4.389996| 0.0003247| 0.0289610| 0.0261347|PMCH/PDYN/SLC6A4/VGF/SYT5/SLC32A1/GAD2/HTR1B/CNR1/CARTPT/SCN2B |    11|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling        |11/46     |495/6858 |  0.0222222|       3.313044| 4.389996| 0.0003247| 0.0289610| 0.0261347|PMCH/PDYN/SLC6A4/VGF/SYT5/SLC32A1/GAD2/HTR1B/CNR1/CARTPT/SCN2B |    11|
## |GO:0010765 |GO:0010765 |positive regulation of sodium ion transport |3/46      |24/6858  |  0.1250000|      18.635870| 7.111704| 0.0005180| 0.0410515| 0.0370452|TESC/SCN3B/SCN2B                                               |     3|
## |GO:0019098 |GO:0019098 |reproductive behavior                       |3/46      |25/6858  |  0.1200000|      17.890435| 6.952064| 0.0005859| 0.0410515| 0.0370452|SLC6A4/PPP1R1B/CNR1                                            |     3|
## |GO:0050433 |GO:0050433 |regulation of catecholamine secretion       |3/46      |25/6858  |  0.1200000|      17.890435| 6.952064| 0.0005859| 0.0410515| 0.0370452|HTR1B/CNR1/CARTPT                                              |     3|
## |GO:0071300 |GO:0071300 |cellular response to retinoic acid          |3/46      |27/6858  |  0.1111111|      16.565217| 6.658917| 0.0007381| 0.0482710| 0.0435602|TNC/SLC6A4/TESC                                                |     3|
## 
## $STR
## 
## 
## |           |ID         |Description                                             |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|    zScore|    pvalue|  p.adjust|    qvalue|geneID                                                                                                             | Count|
## |:----------|:----------|:-------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|---------:|---------:|---------:|:------------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission                          |19/47     |495/6858 |  0.0383838|       5.600774|  8.826999| 0.0000000| 0.0000001| 0.0000001|CHRM1/CAMKV/AKAP5/SLC17A7/GABRG2/KCNIP2/ACHE/NPY/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/SLC1A2/LAMP5 |    19|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling                    |19/47     |495/6858 |  0.0383838|       5.600774|  8.826999| 0.0000000| 0.0000001| 0.0000001|CHRM1/CAMKV/AKAP5/SLC17A7/GABRG2/KCNIP2/ACHE/NPY/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/SLC1A2/LAMP5 |    19|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport         |9/47      |101/6858 |  0.0891089|      13.002317| 10.093909| 0.0000000| 0.0000074| 0.0000060|CHRM1/AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/THY1/SLC9A1                                                           |     9|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission            |14/47     |346/6858 |  0.0404624|       5.904071|  7.775855| 0.0000000| 0.0000094| 0.0000076|CHRM1/CAMKV/AKAP5/ACHE/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/LAMP5                                  |    14|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling                  |14/47     |347/6858 |  0.0403458|       5.887056|  7.760662| 0.0000000| 0.0000094| 0.0000076|CHRM1/CAMKV/AKAP5/ACHE/SYT1/STX1A/GRM2/HOMER1/SV2B/SYNGR1/PLCB1/SYP/PPP1R9A/LAMP5                                  |    14|
## |GO:0051050 |GO:0051050 |positive regulation of transport                        |15/47     |461/6858 |  0.0325380|       4.747773|  6.920638| 0.0000002| 0.0000397| 0.0000321|CHRM1/AKAP5/RAB27B/SCN3B/KCNIP2/ACHE/ACTN2/ATP2B1/SYT1/STX1A/HOMER1/PLCB1/THY1/SLC9A1/SLC1A2                       |    15|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport                  |10/47     |227/6858 |  0.0440529|       6.427969|  6.908302| 0.0000023| 0.0003719| 0.0003007|CHRM1/AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/THY1/SLC9A1                                                     |    10|
## |GO:0010959 |GO:0010959 |regulation of metal ion transport                       |9/47      |188/6858 |  0.0478723|       6.985288|  6.912125| 0.0000040| 0.0005689| 0.0004599|AKAP5/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/THY1/SLC9A1                                                           |     9|
## |GO:0095500 |GO:0095500 |acetylcholine receptor signaling pathway                |4/47      |22/6858  |  0.1818182|      26.529981|  9.962570| 0.0000129| 0.0016365| 0.0013230|CHRM1/ACHE/PLCB1/LY6H                                                                                              |     4|
## |GO:0051928 |GO:0051928 |positive regulation of calcium ion transport            |5/47      |49/6858  |  0.1020408|      14.889275|  8.104878| 0.0000185| 0.0019256| 0.0015567|AKAP5/ATP2B1/HOMER1/THY1/SLC9A1                                                                                    |     5|
## |GO:1905145 |GO:1905145 |cellular response to acetylcholine                      |4/47      |24/6858  |  0.1666667|      24.319149|  9.505860| 0.0000186| 0.0019256| 0.0015567|CHRM1/ACHE/PLCB1/LY6H                                                                                              |     4|
## |GO:0098926 |GO:0098926 |postsynaptic signal transduction                        |4/47      |25/6858  |  0.1600000|      23.346383|  9.297841| 0.0000220| 0.0019300| 0.0015603|CHRM1/ACHE/PLCB1/LY6H                                                                                              |     4|
## |GO:1905144 |GO:1905144 |response to acetylcholine                               |4/47      |25/6858  |  0.1600000|      23.346383|  9.297841| 0.0000220| 0.0019300| 0.0015603|CHRM1/ACHE/PLCB1/LY6H                                                                                              |     4|
## |GO:0098660 |GO:0098660 |inorganic ion transmembrane transport                   |12/47     |458/6858 |  0.0262009|       3.823098|  5.194924| 0.0000414| 0.0033664| 0.0027216|GABRB1/AKAP5/SLC17A7/SCN3B/GABRG2/KCNIP2/ACTN2/ATP2B1/DPP10/PLCB1/THY1/SLC9A1                                      |    12|
## |GO:0034764 |GO:0034764 |positive regulation of transmembrane transport          |6/47      |99/6858  |  0.0606061|       8.843327|  6.529610| 0.0000517| 0.0037049| 0.0029951|AKAP5/KCNIP2/ACTN2/THY1/SLC9A1/SLC1A2                                                                              |     6|
## |GO:0050808 |GO:0050808 |synapse organization                                    |11/47     |394/6858 |  0.0279188|       4.073766|  5.220100| 0.0000521| 0.0037049| 0.0029951|ICAM5/CAMKV/ACTN1/SYNPO/GABRG2/ACHE/ACTN2/HOMER1/PLXNA1/RAP2A/TNC                                                  |    11|
## |GO:1904062 |GO:1904062 |regulation of monoatomic cation transmembrane transport |7/47      |153/6858 |  0.0457516|       6.675845|  5.897764| 0.0000713| 0.0046683| 0.0037740|AKAP5/KCNIP2/ACTN2/HOMER1/DPP10/THY1/SLC9A1                                                                        |     7|
## |GO:1904064 |GO:1904064 |positive regulation of cation transmembrane transport   |5/47      |65/6858  |  0.0769231|      11.224223|  6.879656| 0.0000738| 0.0046683| 0.0037740|AKAP5/KCNIP2/ACTN2/THY1/SLC9A1                                                                                     |     5|
## |GO:0030001 |GO:0030001 |metal ion transport                                     |11/47     |415/6858 |  0.0265060|       3.867624|  5.006252| 0.0000836| 0.0050070| 0.0040478|AKAP5/SLC17A7/SCN3B/KCNIP2/ACTN2/ATP2B1/HOMER1/DPP10/PLCB1/THY1/SLC9A1                                             |    11|
## |GO:1901699 |GO:1901699 |cellular response to nitrogen compound                  |10/47     |351/6858 |  0.0284900|       4.157120|  5.043880| 0.0001027| 0.0058442| 0.0047246|CHRM1/GABRB1/GABRG2/ACHE/ACTN2/ATP2B1/PLCB1/SLC9A1/SLC1A2/LY6H                                                     |    10|
no.sig <- markers.enrich %>% 
  lapply(getElement, "result") %>% 
  sapply(\(x) sum(x$p.adjust <= 0.05))

markers.enrich %>% 
  .[no.sig > 0] %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

DISEASE-SPECIFIC PROTEINS WITH IMPUTATION

Prepare

dat.clean <- dat.clean.orig
dat.clean@colData@listData$area <- dat.clean@colData@listData$condition
dat.clean@colData@listData$areaRep <- dat.clean@colData@listData$replicate
dat.clean@colData@listData$condition <- dat.clean@colData@listData$Diagnosis
dat.clean@colData@listData$replicate <- dat.clean@colData@listData$diagRep
dat.clean@elementMetadata$ID <- dat.clean@elementMetadata$CER_ID
dat.clean@elementMetadata$name <- dat.clean@elementMetadata$CER_name
dat.clean@colData@rownames <- dat.clean$ID

dat.all <- dat.clean

FC

Extract results

min.l2fc = 0.1
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "FC"]

all_contrasts <- test_diff(dat.clean, type = "all")
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))

res <- dep_all@elementMetadata@listData %>%
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

res.split <- comb %>% 
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
                    # filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc) %>% 
                    arrange(!!sym(paste0(cc, "_p.adj")))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1

## 
## $V2

## 
## $V3

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Clustering of DEPs per sample

plot_heatmap(dep_all, type = "centered", kmeans = TRUE, 
             k = 2, col_limit = 3, show_row_names = FALSE, # k 1 or 2?
             indicate = c("condition"))

# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
# 
# rowIds <- ComplexHeatmap::row_order(ht) %>% 
#   lapply(\(x) protIds[x])
# 
# go.res <- rowIds %>% 
#   lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
# 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = TRUE, 
             k = 3, col_limit = 5, show_row_names = FALSE)

Onthology on top significant proteins

markers.enrich <- res.split %>%
  lapply(dplyr::slice, seq(50)) %>% 
  .[sapply(., nrow) > 0] %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

markers.enrich %<>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

markers.enrich %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]

## 
## [[2]]

CER

Extract results

min.l2fc = 0.1
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "CER"]

all_contrasts <- test_diff(dat.clean, type = "all")
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))

res <- dep_all@elementMetadata@listData %>%
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

res.split <- comb %>% 
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
                    # filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc) %>% 
                    arrange(!!sym(paste0(cc, "_p.adj")))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1

## 
## $V2

## 
## $V3

## 
## $V4
## Warning: ggrepel: 28 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Clustering of DEPs per sample

ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE, 
             k = 2, col_limit = 3, show_row_names = FALSE,
             indicate = c("condition"))

protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()

rowIds <- ComplexHeatmap::row_order(ht) %>% 
  lapply(\(x) protIds[x])

go.res <- rowIds %>% 
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

go.res %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
## [[1]]

## 
## [[2]]

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = TRUE, 
             k = 3, col_limit = 5, show_row_names = FALSE)

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: MIURF,XNDC1N,TP53,MRTO4,AURKC,C1QBP
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: RPS28,RPA2,TWNK,MORF4L2,OOEP,KPNA2
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: MEF2A,KPNA1,NBN,RPF2,XRCC1,IL27RA
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: RPSA,MAD2L2,MCAT,HELB,TERF2,MORF4L2
## --> return NULL...
res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: TBX21,VPS72,H1-8,BOP1,LIG4,RAD51
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: KMT5C,IL4,KPNA2,APLF,HELQ,KAT5
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: SMARCAD1,MEAF6,SLC25A36,SLC5A1,SPIDR,STAT6
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: IL4,SPIDR,NEIL3,H1-8,SLC25A33,KDM1A
## --> return NULL...
res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]

## 
## [[2]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]

## 
## [[2]]

STR

Extract results

min.l2fc = 2
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "STR"]

all_contrasts <- test_diff(dat.clean, type = "all")
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))

res <- dep_all@elementMetadata@listData %>%
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

res.split <- comb %>% 
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj")))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1
## Warning: ggrepel: 62 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V2

## 
## $V3

## 
## $V4
## Warning: ggrepel: 37 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5
## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>% 
  ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
  geom_point(size = 3) +
  theme_bw()

Clustering of DEPs per sample

ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE, 
             k = 2, col_limit = 3, show_row_names = F,
             indicate = c("condition"))

protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()

rowIds <- ComplexHeatmap::row_order(ht) %>% 
  lapply(\(x) protIds[x])

go.res <- rowIds %>% 
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

go.res %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
## [[1]]

## 
## [[2]]

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = TRUE, 
             k = 3, col_limit = 5, show_row_names = FALSE)

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: HDGFL2,MIURF,XNDC1N,TP53,MRTO4,AURKC
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: C1QBP,RPS28,RPA2,TWNK,MORF4L2,OOEP
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: KPNA2,MEF2A,KPNA1,NBN,RPF2,XRCC1
## --> return NULL...
res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: IL27RA,RPSA,MAD2L2,MCAT,HELB,TERF2
## --> return NULL...
res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]

## 
## [[2]]

STRINGdb network

Networks are shown for all evidence of interaction, not just physical evidence (stronger).

sdb <- STRINGdb::STRINGdb(species = 9606, score_threshold = 400, version = "12.0")
sdb.prots <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    message(x)
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_"))) %>%
      pull(name)
  }) %>% 
  setNames(names(res.split)) %>%
  .[sapply(., length) > 1]
## PSP_vs_MSA
## PSP_vs_PD
## PSP_vs_CTRL
## MSA_vs_PD
## MSA_vs_CTRL
## PD_vs_CTRL
sdb.prots %>% 
  lapply(sdb$plot_network)

## $PSP_vs_MSA
## NULL
## 
## $MSA_vs_PD
## NULL
## 
## $MSA_vs_CTRL
## NULL

SN

Extract results

min.l2fc = 0.5
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "SN"]

all_contrasts <- test_diff(dat.clean, type = "all")
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))

res <- dep_all@elementMetadata@listData %>%
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

res.split <- comb %>% 
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj")))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)
## $V1

## 
## $V2

## 
## $V3
## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 100, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Clustering of DEPs per sample

ht <- plot_heatmap(dep_all, type = "centered", kmeans = TRUE, 
             k = 2, col_limit = 3, show_row_names = FALSE,
             indicate = c("condition"))

protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()

rowIds <- ComplexHeatmap::row_order(ht) %>% 
  lapply(\(x) protIds[x])

go.res <- rowIds %>% 
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

go.res %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))
## [[1]]

## 
## [[2]]

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = TRUE, 
             k = 2, col_limit = 5, show_row_names = FALSE)

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: CGAS,KMT5B,HDGFL2,MIURF,XNDC1N,TP53
## --> return NULL...
res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]
## --> No gene can be mapped....
## --> Expected input gene ID: MRTO4,AURKC,C1QBP,RPS28,RPA2,TWNK
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: MORF4L2,OOEP,KPNA2,MEF2A,KPNA1,NBN
## --> return NULL...
## --> No gene can be mapped....
## --> Expected input gene ID: RPF2,XRCC1,IL27RA,RPSA,MAD2L2,MCAT
## --> return NULL...
res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))
## [[1]]

## 
## [[2]]

MED

Extract results

min.l2fc = 0.1
min.padj = 0.2

dat.clean <- dat.all[, dat.all$area == "MED"]

all_contrasts <- test_diff(dat.clean, type = "all")
## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL
comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))

res <- dep_all@elementMetadata@listData %>%
  as.data.frame() %>% 
  mutate(imputed = apply(.[, grepl("imputed", colnames(.))], 1, \(x) any(x)),
         num_NAs = apply(.[, grepl("num_NAs", colnames(.))], 1, \(x) sum(x))) %>%
  filter(num_NAs < 121) # 70 % = 121, 50 % = 87, 25 % = 43

res.split <- comb %>% 
  setNames(lapply(., \(cc) res[, c("name", "imputed", "num_NAs", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj")))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)
## $V1

## 
## $V2

## 
## $V3

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")
## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.
## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Onthology on top significant proteins

markers.enrich <- res.split %>%
  lapply(dplyr::slice, seq(50)) %>% 
  .[sapply(., nrow) > 0] %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T)

markers.enrich %<>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

markers.enrich %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Session

Sys.time() - tt
## Time difference of 11.93094 mins
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2;  LAPACK version 3.12.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DEP_1.26.0                  qs_0.27.3                  
##  [3] SummarizedExperiment_1.38.1 Biobase_2.68.0             
##  [5] GenomicRanges_1.60.0        GenomeInfoDb_1.44.0        
##  [7] IRanges_2.42.0              S4Vectors_0.46.0           
##  [9] BiocGenerics_0.54.0         generics_0.1.4             
## [11] MatrixGenerics_1.20.0       matrixStats_1.5.0          
## [13] pheatmap_1.0.12             biomaRt_2.64.0             
## [15] scHelper_0.0.5              ggpubr_0.6.0               
## [17] ggrepel_0.9.6               cowplot_1.1.3              
## [19] ggplot2_3.5.2               dplyr_1.1.4                
## [21] magrittr_2.0.3             
## 
## loaded via a namespace (and not attached):
##   [1] STRINGdb_2.20.0             bitops_1.0-9               
##   [3] fs_1.6.6                    ProtGenerics_1.40.0        
##   [5] enrichplot_1.28.2           httr_1.4.7                 
##   [7] RColorBrewer_1.1-3          doParallel_1.0.17          
##   [9] ggsci_3.2.0                 tools_4.5.2                
##  [11] MSnbase_2.30.1              backports_1.5.0            
##  [13] R6_2.6.1                    DT_0.33                    
##  [15] lazyeval_0.2.2              GetoptLong_1.0.5           
##  [17] withr_3.0.2                 prettyunits_1.2.0          
##  [19] gridExtra_2.3               preprocessCore_1.70.0      
##  [21] fdrtool_1.2.18              cli_3.6.5                  
##  [23] Cairo_1.6-2                 sandwich_3.1-1             
##  [25] labeling_0.4.3              conos_1.5.2                
##  [27] sass_0.4.10                 mvtnorm_1.3-3              
##  [29] readr_2.1.5                 yulab.utils_0.2.0          
##  [31] gson_0.1.0                  DOSE_4.2.0                 
##  [33] R.utils_2.13.0              dichromat_2.0-0.1          
##  [35] MetaboCoreUtils_1.16.1      plotrix_3.8-4              
##  [37] limma_3.64.0                rstudioapi_0.17.1          
##  [39] impute_1.82.0               RSQLite_2.3.11             
##  [41] gridGraphics_0.5-1          shape_1.4.6.1              
##  [43] RApiSerialize_0.1.4         gtools_3.9.5               
##  [45] car_3.1-3                   GO.db_3.21.0               
##  [47] Matrix_1.7-3                MALDIquant_1.22.3          
##  [49] imputeLCMD_2.1              abind_1.4-8                
##  [51] R.methodsS3_1.8.2           lifecycle_1.0.4            
##  [53] yaml_2.3.10                 carData_3.0-5              
##  [55] gplots_3.2.0                qvalue_2.40.0              
##  [57] SparseArray_1.8.0           BiocFileCache_2.16.0       
##  [59] Rtsne_0.17                  grid_4.5.2                 
##  [61] blob_1.2.4                  promises_1.3.2             
##  [63] crayon_1.5.3                PSMatch_1.12.0             
##  [65] shinydashboard_0.7.3        ggtangle_0.0.6             
##  [67] lattice_0.22-7              mzR_2.38.0                 
##  [69] KEGGREST_1.48.0             magick_2.8.6               
##  [71] pillar_1.10.2               knitr_1.50                 
##  [73] ComplexHeatmap_2.24.0       fgsea_1.34.0               
##  [75] rjson_0.2.23                codetools_0.2-20           
##  [77] fastmatch_1.1-6             glue_1.8.0                 
##  [79] ggfun_0.1.8                 pcaMethods_2.0.0           
##  [81] data.table_1.17.2           MultiAssayExperiment_1.34.0
##  [83] vctrs_0.6.5                 png_0.1-8                  
##  [85] treeio_1.32.0               gtable_0.3.6               
##  [87] gsubfn_0.7                  assertthat_0.2.1           
##  [89] cachem_1.1.0                xfun_0.52                  
##  [91] S4Arrays_1.8.0              mime_0.13                  
##  [93] ncdf4_1.24                  iterators_1.0.14           
##  [95] gmm_1.8                     statmod_1.5.0              
##  [97] nlme_3.1-168                ggtree_3.16.0              
##  [99] bit64_4.6.0-1               progress_1.2.3             
## [101] filelock_1.0.3              bslib_0.9.0                
## [103] affyio_1.78.0               tmvtnorm_1.6               
## [105] KernSmooth_2.23-26          colorspace_2.1-1           
## [107] DBI_1.2.3                   tidyselect_1.2.1           
## [109] chron_2.3-62                bit_4.6.0                  
## [111] compiler_4.5.2              curl_7.0.0                 
## [113] httr2_1.2.1                 xml2_1.4.0                 
## [115] DelayedArray_0.34.1         stringfish_0.16.0          
## [117] caTools_1.18.3              scales_1.4.0               
## [119] affy_1.86.0                 rappdirs_0.3.3             
## [121] stringr_1.5.1               digest_0.6.37              
## [123] rmarkdown_2.29              XVector_0.48.0             
## [125] htmltools_0.5.8.1           pkgconfig_2.0.3            
## [127] dbplyr_2.5.1                fastmap_1.2.0              
## [129] rlang_1.1.6                 GlobalOptions_0.1.2        
## [131] htmlwidgets_1.6.4           UCSC.utils_1.4.0           
## [133] shiny_1.10.0                farver_2.1.2               
## [135] jquerylib_0.1.4             zoo_1.8-14                 
## [137] jsonlite_2.0.0              BiocParallel_1.42.0        
## [139] mzID_1.46.0                 GOSemSim_2.34.0            
## [141] R.oo_1.27.1                 Formula_1.2-5              
## [143] GenomeInfoDbData_1.2.14     ggplotify_0.1.2            
## [145] patchwork_1.3.0             Rcpp_1.0.14                
## [147] proto_1.0.0                 ape_5.8-1                  
## [149] MsCoreUtils_1.20.0          sqldf_0.4-11               
## [151] vsn_3.76.0                  leidenAlg_1.1.5            
## [153] stringi_1.8.7               MASS_7.3-65                
## [155] org.Hs.eg.db_3.21.0         plyr_1.8.9                 
## [157] parallel_4.5.2              Biostrings_2.76.0          
## [159] sccore_1.0.5                splines_4.5.2              
## [161] hash_2.2.6.3                hms_1.1.3                  
## [163] Spectra_1.18.0              circlize_0.4.16            
## [165] igraph_2.1.4                QFeatures_1.18.0           
## [167] ggsignif_0.6.4              reshape2_1.4.4             
## [169] XML_3.99-0.18               evaluate_1.0.3             
## [171] RcppParallel_5.1.10         BiocManager_1.30.25        
## [173] tzdb_0.5.0                  foreach_1.5.2              
## [175] httpuv_1.6.16               tidyr_1.3.1                
## [177] purrr_1.0.4                 clue_0.3-66                
## [179] norm_1.0-11.1               BiocBaseUtils_1.10.0       
## [181] broom_1.0.8                 xtable_1.8-4               
## [183] AnnotationFilter_1.32.0     tidytree_0.4.6             
## [185] rstatix_0.7.2               later_1.4.2                
## [187] tibble_3.2.1                clusterProfiler_4.16.0     
## [189] aplot_0.2.5                 memoise_2.0.1              
## [191] AnnotationDbi_1.70.0        cluster_2.1.8.1